Mapping potato late blight from UAV-based multispectral imagery

Part 1

Jorge Luis Rodríguez

email: jorodriguezga@unal.edu.co

Introduction

This work proposes a method to process high resolution multispectral images for detection and mapping late blight in potato crops. It is based on three main aspects: (i) radiometric and geometric correction of the raw images; (ii) morphological operations for background removing; (iii) model trainig and supervised classification stage using a random forest classifier; and (iv) implementation of trained model in a new dataset.

Data and methods

The study area is a 1920 sq.m. potato field located in Subachoque, Colombia (see figure). This field is part of an experimental plot designed to evaluate the potato Diacol capiro variety response to different nutrient treatments. As the weather conditions favoured the appearance of late blight disease in several plots, the project had the chance to monitor the crop disease development. The experimental plot was inspected every seven days during the entire crop life span (i.e. 120 days from planting to maturity).

Mosaic generated extend over 3.2 Ha, which include the experimental plot and a bigger area where a different variety of potato was located, being this variety a yellow potato, this means, the tuber root colour it produces is yellow. In contrast, the potato from the experimental plot was from one variety of white potato. The experimental plot area was clipped by using a reference polygon (red polygon).

Study area:

Imagen

High resolution multispectral images were acquired at 40 m altitude above the ground surface at 11:00 am local time (GMT-5). Each multispectral image acquired by the MicaSense camera had five bands as described in the next table. In this work camera bands 4 and 5 have been reset according to the following order: Blue (B); Green (G); Red (R); Red edge (RE) and Near infrared (NIR).

# Banda Nombre Centro de Banda (nm) Ancho Banda FWHM (nm)
1 Azul 475 20
2 Verde 560 20
3 Rojo 668 10
4 Infrarrojo cercano
840 40
5 Borde Rojo 717 10

Experimental crop area consisted of an array of 18 blocks of 12 m $\times$ 8 m inside of a field of 77 m $\times$ 24 m. Each block had 9 rows arranged along the long side of the field with an area of 216 sq.m., the space between rows was 1m. Each row contained approximately 30 seed tubers of the same variety.

Imagen (a) May 12, 2018 Imagen (b) May 26, 2018

The layout of the experimental crop. The yellow lines indicate the division of the field into 18 experimental blocks. (a) May 12, 2018; (b) May 26, 2018.

Background removing

Multispectral image load

In [1]:
%matplotlib inline
import time
import os, shutil
from matplotlib import pyplot as plt
from matplotlib import colors
from IPython.display import Image
import numpy as np
import cv2
import math
import time
from osgeo import gdal
from skimage import img_as_ubyte
import scipy
import pandas as pd
from skimage.filters import threshold_otsu
import skimage.io as io
  
from sklearn.ensemble import AdaBoostClassifier, RandomForestClassifier, GradientBoostingClassifier, ExtraTreesClassifier
from sklearn.externals import joblib
from sklearn.metrics import plot_roc_curve, classification_report, multilabel_confusion_matrix
from sklearn.model_selection import train_test_split

from sklearn.svm import SVC, LinearSVC
from sklearn.neighbors import KNeighborsClassifier
/usr/local/lib/python3.5/dist-packages/sklearn/externals/joblib/__init__.py:15: FutureWarning: sklearn.externals.joblib is deprecated in 0.21 and will be removed in 0.23. Please import this functionality directly from joblib, which can be installed with: pip install joblib. If this warning is raised when loading pickled models, you may need to re-serialize those models with scikit-learn 0.21+.
  warnings.warn(msg, category=FutureWarning)
In [2]:
ds = gdal.Open('Raster/Correccion_reflectancia/MR_CE_20180512_Subachoque.tif')
In [3]:
# asignacion cada banda 
AZUL = ds.GetRasterBand(1).ReadAsArray()
VERDE = ds.GetRasterBand(2).ReadAsArray()
ROJO = ds.GetRasterBand(3).ReadAsArray()
REDEDGE = ds.GetRasterBand(4).ReadAsArray()
NIR = ds.GetRasterBand(5).ReadAsArray()

srs = ds.GetProjectionRef()
geo_transform = ds.GetGeoTransform()

plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(321) ,plt.imshow(AZUL, cmap='gray'),plt.title('Blue band')
plt.subplot(322) ,plt.imshow(VERDE, cmap='gray'),plt.title('Green band')
plt.subplot(323) ,plt.imshow(ROJO, cmap='gray'),plt.title('Red band')
plt.subplot(324) ,plt.imshow(REDEDGE, cmap='gray'),plt.title('Red edge band')
plt.subplot(325) ,plt.imshow(NIR, cmap='gray'),plt.title('Infrared band')
plt.show()

Thresholding

A thresholding technique reduces a grey-level image into an image where objects and background are represented with two levels: a binary image (Glasbey, 1993). Because of the difference in reflectance between the soil and the potato plants, both near infrared and red edge bands allowed separation of the potato plants from the soil.

The first step was to analyse multispectral bands of the orthomosaics by plotting their histograms. The next figure shows the histogram for each band in the multispectral orthomosaic for the two categories under study as well as ground reference images of the two categories. Green, Rededge and Near infrared bands have right-skewed bimodal histograms in the first category. Bands with high weed presence had a right-skewed histograms without a clear valley. Thresholding method used to separate vegetation was OTSU's thresholding method (Otsu, 1979).

In [4]:
hist_b1 = np.histogram(AZUL[AZUL!=0].ravel(),100,[0,1])[0]
hist_b2 = np.histogram(VERDE[VERDE!=0].ravel(),100,[0,1])[0]
hist_b3 = np.histogram(ROJO[ROJO!=0].ravel(),100,[0,1])[0]
hist_b4 = np.histogram(REDEDGE[REDEDGE!=0].ravel(),100,[0,1])[0]
hist_b5 = np.histogram(NIR[NIR!=0].ravel(),100,[0,1])[0]

hist_range = np.arange(0,1,0.01)
plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(231), plt.bar(hist_range,hist_b1,0.01), plt.title("Blue")
plt.subplot(232), plt.bar(hist_range,hist_b2,0.01),plt.title('Green')
plt.subplot(233), plt.bar(hist_range,hist_b3,0.01),plt.title('Red')
plt.subplot(234), plt.bar(hist_range,hist_b4,0.01),plt.title('Rededge')
plt.subplot(235), plt.bar(hist_range,hist_b5,0.01),plt.title('Near infrared')
Out[4]:
(<matplotlib.axes._subplots.AxesSubplot at 0x7f99b2446390>,
 <Container object of 100 artists>,
 <matplotlib.text.Text at 0x7f99bdcb04a8>)
In [5]:
image=NIR
hist_image = image.ravel()
hist = np.histogram(image[image!=0].ravel(),100,[0,1])[0]
thresh = threshold_otsu(image)
binary = image > thresh

plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(131) ,plt.imshow(NIR, cmap='gray'),plt.title('Near infrared')
plt.subplot(132) ,plt.imshow(binary,cmap='gray'),plt.title('Binary image')
plt.show()

plt.show()

Masking image

To delete background from original multispectral image we multiply each band for the binary image resulting from thresholding step using the next equation:

\begin{equation} \label{Eq3} C_{i}(x,y) = A_{i}(x,y)B(x,y), \end{equation}

where $C_{i}(x,y)$ is the resulting pixel at position $(x,y)$ for $i$ band without background, $A_{i}(x,y)$ is the original pixel at position $(x,y)$ for $i$ band of the multispectral image and $B(x,y)$ is the binary image pixel at position $(x,y)$ created in the thresholding process.

In [6]:
mask = np.array(binary, dtype='uint8')
filtered_mask = cv2.medianBlur(mask,3)

plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(131) ,plt.imshow(NIR, cmap='gray'),plt.title('Banda NIR')
plt.subplot(132) ,plt.imshow(mask,cmap='gray'),plt.title('Mascara')
plt.subplot(133) ,plt.imshow(filtered_mask,cmap='gray'),plt.title('Mascara con filtro de mediana')
plt.show()

plt.show()

Here we can see how looks any band from the multispectral image after masking process. It can be seen how the ground surface does not appear any more.

In [7]:
azul_seg = AZUL*filtered_mask
verde_seg = VERDE*filtered_mask
rojo_seg = ROJO*filtered_mask
rededge_seg = REDEDGE*filtered_mask
nir_seg = NIR*filtered_mask

plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(321) ,plt.imshow(azul_seg, cmap='gray'),plt.title('Blue')
plt.subplot(322) ,plt.imshow(verde_seg, cmap='gray'),plt.title('Green')
plt.subplot(323) ,plt.imshow(rojo_seg, cmap='gray'),plt.title('Red')
plt.subplot(324) ,plt.imshow(rededge_seg, cmap='gray'),plt.title('Red edge')
plt.subplot(325) ,plt.imshow(nir_seg, cmap='gray'),plt.title('Near infrared')
plt.show()

We can now create a false colour composition to improve the visualization of the background removing step. This was performed for visualisation purposes only:

In [8]:
# Contrast stretch function for enhance the bands
# Only for visualization purposes
def enhance_band(band, min_in= 0.1, max_in= 0.2):
    min_setup = 0.0
    max_setup = 1.0
    enhanced = (band-min_in)*(((max_setup-min_setup)/(max_in-min_in))+min_setup)
    return enhanced
In [9]:
NIR_IN = enhance_band(NIR, min_in= 0.0, max_in= (NIR.max()-0.12))
ROJO_IN = enhance_band(ROJO, min_in= 0.0, max_in= (ROJO.max()-0.06))
VERDE_IN = enhance_band(VERDE, min_in= 0.0, max_in= VERDE.max())

rgb = np.stack([NIR_IN,ROJO_IN,VERDE_IN], axis=2) #axis =2 so its shape is (M,N,3) otherwise doesn't work with plt. (M, N, 3): an image with RGB values (0-1 float or 0-255 int)
mask = cv2.cvtColor(filtered_mask, cv2.COLOR_GRAY2BGR) # So we can apply the mask 
plantsImage = rgb*mask

plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(121) ,plt.imshow(rgb),plt.title("Original image")
plt.subplot(122) ,plt.imshow(plantsImage),plt.title("Masked image")
plt.show()

plt.figure(2, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(121) ,plt.imshow(rgb[1450:1700,1100:1300]),plt.title("Original image")
plt.subplot(122) ,plt.imshow(plantsImage[1450:1700,1100:1300]),plt.title("Masked image")
plt.show()

Segmented bands to files

Now we proceed to save each of the segmented bands:

In [10]:
bands_out = np.stack([azul_seg,
                      verde_seg,
                      rojo_seg,
                      rededge_seg,
                     nir_seg],
                     axis=2)

band_names = {'M_2018-05-12_Corte_Seg_B1.tif':0,
              'M_2018-05-12_Corte_Seg_B2.tif':1,
             'M_2018-05-12_Corte_Seg_B3.tif':2,
             'M_2018-05-12_Corte_Seg_B4.tif':3,
             'M_2018-05-12_Corte_Seg_B5.tif':4
             }
dir_out = 'Raster/Segmented/'

for band_name in band_names:
    filename_output = dir_out + band_name
    index = band_names[band_name]
    x_size = ds.RasterXSize  # Raster xsize
    y_size = ds.RasterYSize  # Raster ysize
    driver = gdal.GetDriverByName('GTiff')
    arch = driver.Create(filename_output,x_size,y_size,1,gdal.GDT_Float32)
    arch.SetGeoTransform(geo_transform)
    arch.SetProjection(srs)
    arch.GetRasterBand(1).WriteArray(bands_out[:,:,index].astype(np.float32))
    del(arch)
    print("Band "+band_name+" exported")
Band M_2018-05-12_Corte_Seg_B5.tif exported
Band M_2018-05-12_Corte_Seg_B4.tif exported
Band M_2018-05-12_Corte_Seg_B1.tif exported
Band M_2018-05-12_Corte_Seg_B3.tif exported
Band M_2018-05-12_Corte_Seg_B2.tif exported

Supervised classification of dataset (a) May 12, 2018, and models training

Loading segmented bands

In [11]:
rootdir = "Raster/Segmented"
# path to your training data
path_pix = "Entrenamiento/"
# path to your model
path_model = "Data/Models/"
# path to your classification results
path_class = "Data/Class/"

samples = path_pix + "Entrenamiento_raster.tif"

ds_seg1 = gdal.Open('Raster/Segmented/M_2018-05-12_Corte_Seg_B1.tif')
ds_seg2 = gdal.Open('Raster/Segmented/M_2018-05-12_Corte_Seg_B2.tif')
ds_seg3 = gdal.Open('Raster/Segmented/M_2018-05-12_Corte_Seg_B3.tif')
ds_seg4 = gdal.Open('Raster/Segmented/M_2018-05-12_Corte_Seg_B4.tif')
ds_seg5 = gdal.Open('Raster/Segmented/M_2018-05-12_Corte_Seg_B5.tif')

Segmented_image = np.stack([ds_seg1.GetRasterBand(1).ReadAsArray(),
                           ds_seg2.GetRasterBand(1).ReadAsArray(),
                           ds_seg3.GetRasterBand(1).ReadAsArray(),
                           ds_seg4.GetRasterBand(1).ReadAsArray(),
                           ds_seg5.GetRasterBand(1).ReadAsArray()], axis=2)

Training and classification of the dataset is performed by the next function:

In [12]:
def plot_time(dt):
    fig, ax = plt.subplots()

    methods = ('RF', 'GBC', 'SVC', 'LSVC', 'KNN')
    y_pos = np.arange(len(methods))
    performance = dt

    ax.barh(y_pos, performance, align='center')
    ax.set_yticks(y_pos)
    ax.set_yticklabels(methods)
    ax.invert_yaxis()  # labels read top-to-bottom
    ax.set_xlabel('Seconds')
    ax.set_title('Performance by method')

    plt.show()

# declare a new function
def training(img_ds, samples, path_pix, path_model, path_class):

    img = img_ds.copy()

    # Load trainig classes 
    roi_ds = io.imread(samples)   
    roi = np.array(roi_ds, dtype='int8')  
    
    # read in your labels
    labels = np.unique(roi[roi > 0]) 
    print('The training data include {n} classes: {classes}'.format(n=labels.size, classes=labels))

    # compose your X,Y data (dataset - training data)     
    X = img[roi > 0, :] 
    Y = roi[roi > 0] 
    
    # Split dataset into trainig an test
    X_train, X_test, Y_train, Y_test = train_test_split(X, Y, random_state=42)

    # assign class weights (class 1 has the weight 3, etc.)
    #weights = {1:3, 2:2, 3:2, 4:2}
    weights = None
    #weights = "balanced_subsample"
    
    ################### Random Forest Classifier ##########################
    # build your Random Forest Classifier 
    # for more information: http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html
    ti_rf = time.time()
    rf = RandomForestClassifier(class_weight = weights, n_estimators = 500, criterion = 'gini', max_depth = 4, 
                                min_samples_split = 2, min_samples_leaf = 2, max_features = 'auto', 
                                bootstrap = True, oob_score = True, n_jobs = 2, random_state = None, verbose = True)  
    
    
    # now fit your training data with the original dataset
    rf = rf.fit(X_train,Y_train)

    # Save your Random Forest Model     
    model = path_model + "model_RF.pkl"
    joblib.dump(rf, model)   

    # call your random forest model
    rf = path_model + "model_RF.pkl"          
    clf = joblib.load(rf)    

    # Classification of array and save as image (23 refers to the number of multitemporal NDVI bands in the stack) 
    new_shape = (img.shape[0] * img.shape[1], img.shape[2]) 
    img_as_array = img[:, :, :23].reshape(new_shape)   

    class_prediction = clf.predict(img_as_array) 
    class_prediction = class_prediction.reshape(img[:, :, 0].shape)     
    tf_rf = time.time()
    dt_rf = tf_rf - ti_rf
    
    #########################Gradient Boost Classifier#############################
    # alternatively you may try out a Gradient Boosting Classifier 
    # It is much less RAM consuming and considers weak training data      
    ti_gbc = time.time()
    GBC = GradientBoostingClassifier(n_estimators = 500, min_samples_leaf = 1, min_samples_split = 4, max_depth = 4,    
                                    max_features = 'auto', learning_rate = 0.8, subsample = 1, random_state = None,         
                                    warm_start = True)
    
    # now fit your training data with the original dataset
    GBC = GBC.fit(X_train,Y_train)

    # Save your Gradient Boosting Model     
    model_GBC = path_model + "model_GBC.pkl"
    joblib.dump(GBC, model_GBC)

    GBC = path_model + "model_GBC.pkl"          
    clf_GBC = joblib.load(GBC)    

    # Classification of array and save as image (23 refers to the number of multitemporal NDVI bands in the stack) 
    new_shape = (img.shape[0] * img.shape[1], img.shape[2]) 
    img_as_array = img[:, :, :23].reshape(new_shape)   

    class_prediction_GBC = clf_GBC.predict(img_as_array) 
    class_prediction_GBC = class_prediction_GBC.reshape(img[:, :, 0].shape)
    tf_gbc = time.time()
    dt_gbc = tf_gbc - ti_gbc
        
    #########################Support Vector Classifier#############################    
    ti_svc = time.time()
    svc = SVC(random_state = None, gamma = 'auto')
    
    # now fit your training data with the original dataset
    svc = svc.fit(X_train,Y_train)

    # export your SVC model     
    model_svc = path_model + "model_svc.pkl"
    joblib.dump(svc, model_svc)

    # call your SVC model
    svc = path_model + "model_svc.pkl"          
    clf_svc = joblib.load(svc)    

    # Classification of array and save as image (23 refers to the number of multitemporal NDVI bands in the stack) 
    new_shape = (img.shape[0] * img.shape[1], img.shape[2]) 
    img_as_array = img[:, :, :23].reshape(new_shape)   

    class_prediction_svc = clf_svc.predict(img_as_array) 
    class_prediction_svc = class_prediction_svc.reshape(img[:, :, 0].shape)  
    tf_svc = time.time()
    dt_svc = tf_svc - ti_svc
    #########################Linear Support Vector Classifier#############################   
    ti_lsvc = time.time()
    lsvc = LinearSVC()
    
    # now fit your training data with the original dataset
    lsvc = lsvc.fit(X_train,Y_train)

    # export your Random Forest / Gradient Boosting Model     
    model_lsvc = path_model + "model_lsvc.pkl"
    joblib.dump(lsvc, model_lsvc)

    # call your Linear SVC model
    lsvc = path_model + "model_lsvc.pkl"          
    clf_lsvc = joblib.load(lsvc)    

    # Classification of array and save as image (23 refers to the number of multitemporal NDVI bands in the stack) 
    new_shape = (img.shape[0] * img.shape[1], img.shape[2]) 
    img_as_array = img[:, :, :23].reshape(new_shape)   

    class_prediction_lsvc = clf_lsvc.predict(img_as_array) 
    class_prediction_lsvc = class_prediction_lsvc.reshape(img[:, :, 0].shape)  

    tf_lsvc = time.time()
    dt_lsvc = tf_lsvc - ti_lsvc
    ######################### KNN Classifier#############################  
    ti_knn = time.time()
    knn = KNeighborsClassifier()
    
    # now fit your training data with the original dataset
    knn = knn.fit(X_train,Y_train)

    # export your KNN Model     
    model_knn = path_model + "model_knn.pkl"
    joblib.dump(knn, model_knn)

    # call your KNN model
    knn = path_model + "model_knn.pkl"          
    clf_knn = joblib.load(knn)    

    # Classification of array and save as image (23 refers to the number of multitemporal NDVI bands in the stack) 
    new_shape = (img.shape[0] * img.shape[1], img.shape[2]) 
    img_as_array = img[:, :, :23].reshape(new_shape)   

    class_prediction_knn = clf_knn.predict(img_as_array) 
    class_prediction_knn = class_prediction_knn.reshape(img[:, :, 0].shape) 
    
    tf_knn = time.time()
    dt_knn = tf_knn - ti_knn
    ########## Plot ROC curves ####################
    ax = plt.gca()

    rfc_disp = plot_roc_curve(clf,X_test,Y_test, ax=ax)
    GBC_disp = plot_roc_curve(clf_GBC,X_test,Y_test, ax=rfc_disp.ax_)
    svc_disp = plot_roc_curve(clf_svc,X_test,Y_test, ax=GBC_disp.ax_, alpha=0.8)
    lsvc_disp = plot_roc_curve(clf_lsvc,X_test,Y_test, ax=svc_disp.ax_, alpha=0.8)
    knn_disp = plot_roc_curve(clf_knn,X_test,Y_test, ax=lsvc_disp.ax_, alpha=0.8)
    knn_disp.figure_.suptitle("ROC curve comparison")
    plt.show()
    
    
    ax1 = plt.gca()
    ax1.set(xlim=[0,0.7], ylim=[0.4,1])

    rfc_disp = plot_roc_curve(clf,X_test,Y_test, ax=ax1)
    GBC_disp = plot_roc_curve(clf_GBC,X_test,Y_test, ax=rfc_disp.ax_)
    svc_disp = plot_roc_curve(clf_svc,X_test,Y_test, ax=GBC_disp.ax_, alpha=0.8)
    lsvc_disp = plot_roc_curve(clf_lsvc,X_test,Y_test, ax=svc_disp.ax_, alpha=0.8)
    knn_disp = plot_roc_curve(clf_knn,X_test,Y_test, ax=lsvc_disp.ax_, alpha=0.8)
    knn_disp.figure_.suptitle("ROC curve comparison - ZOOM")
    plt.show()
    
    plot_time([dt_rf, dt_gbc, dt_svc, dt_lsvc,dt_knn])
    
    
    return class_prediction, class_prediction_GBC, class_prediction_svc, class_prediction_lsvc, class_prediction_knn
In [13]:
class_prediction, class_prediction_GBC, class_prediction_svc, class_prediction_lsvc, class_prediction_knn = training(Segmented_image, samples, path_pix, path_model, path_class)
The training data include 2 classes: [1 2]
[Parallel(n_jobs=2)]: Using backend ThreadingBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done  46 tasks      | elapsed:    0.2s
[Parallel(n_jobs=2)]: Done 196 tasks      | elapsed:    0.7s
[Parallel(n_jobs=2)]: Done 446 tasks      | elapsed:    1.5s
[Parallel(n_jobs=2)]: Done 500 out of 500 | elapsed:    1.8s finished
[Parallel(n_jobs=2)]: Using backend ThreadingBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done  46 tasks      | elapsed:    4.2s
[Parallel(n_jobs=2)]: Done 196 tasks      | elapsed:   17.5s
[Parallel(n_jobs=2)]: Done 446 tasks      | elapsed:   38.2s
[Parallel(n_jobs=2)]: Done 500 out of 500 | elapsed:   42.6s finished
[Parallel(n_jobs=2)]: Using backend ThreadingBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done  46 tasks      | elapsed:    0.0s
[Parallel(n_jobs=2)]: Done 196 tasks      | elapsed:    0.0s
[Parallel(n_jobs=2)]: Done 446 tasks      | elapsed:    0.1s
[Parallel(n_jobs=2)]: Done 500 out of 500 | elapsed:    0.1s finished
[Parallel(n_jobs=2)]: Using backend ThreadingBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done  46 tasks      | elapsed:    0.0s
[Parallel(n_jobs=2)]: Done 196 tasks      | elapsed:    0.0s
[Parallel(n_jobs=2)]: Done 446 tasks      | elapsed:    0.1s
[Parallel(n_jobs=2)]: Done 500 out of 500 | elapsed:    0.1s finished

The scikit-klearn methods do not support NaN data, so we can't replace pixels with zero values by using np.nan, so when classifying, pixels with a zero value are included in one of the classes, in this case, in the late blight class. To clean the results, we use the mask created previously.

In [14]:
from matplotlib.colors import from_levels_and_colors
cmap, norm = from_levels_and_colors([0,0.5,1,2],['black','red','yellow'])

plt.figure(3, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(231) ,plt.imshow(class_prediction, cmap=cmap),plt.title('RFC')
plt.subplot(232) ,plt.imshow(class_prediction_GBC,cmap=cmap),plt.title('GBC')
plt.subplot(233) ,plt.imshow(class_prediction_svc,cmap=cmap),plt.title('SVC')
plt.subplot(234) ,plt.imshow(class_prediction_lsvc,cmap=cmap),plt.title('LSVC')
plt.subplot(235) ,plt.imshow(class_prediction_knn,cmap=cmap),plt.title('KNN')
plt.show()
In [15]:
cmap, norm = from_levels_and_colors([0,0.5,1,2],['white','#EA8A00','green'])

rf_classification = class_prediction*filtered_mask
gbc_classification = class_prediction_GBC*filtered_mask
svc_classification = class_prediction_svc*filtered_mask
lsvc_classification = class_prediction_lsvc*filtered_mask
knn_classification = class_prediction_knn*filtered_mask

plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(231) ,plt.imshow(rf_classification, cmap=cmap),plt.title('RFC')
plt.subplot(232) ,plt.imshow(gbc_classification,cmap=cmap),plt.title('GBC')
plt.subplot(233) ,plt.imshow(svc_classification,cmap=cmap),plt.title('SVC')
plt.subplot(234) ,plt.imshow(lsvc_classification,cmap=cmap),plt.title('LSVC')
plt.subplot(235) ,plt.imshow(knn_classification,cmap=cmap),plt.title('KNN')
plt.show()

plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(231) ,plt.imshow(rf_classification[1400:1700,1200:1400], cmap=cmap),plt.title('RFC Zoom')
plt.subplot(232) ,plt.imshow(gbc_classification[1400:1700,1200:1400],cmap=cmap),plt.title('GBC  Zoom')
plt.subplot(233) ,plt.imshow(svc_classification[1400:1700,1200:1400],cmap=cmap),plt.title('SVC  Zoom')
plt.subplot(234) ,plt.imshow(lsvc_classification[1400:1700,1200:1400],cmap=cmap),plt.title('LSVC  Zoom')
plt.subplot(235) ,plt.imshow(knn_classification[1400:1700,1200:1400],cmap=cmap),plt.title('KNN  Zoom')
plt.show()
In [16]:
labels_class = np.unique(class_prediction[class_prediction > 0]) 
print('The classified data include {n} classes: {classes}'.format(n=labels_class.size, classes=labels_class))
The classified data include 2 classes: [1 2]

Accuracy assessment of results

Call ground reference

In [17]:
ds_ground_truth = gdal.Open('Validacion/Raster/Ground_truth_2018-05-12_.tif')

# asignacion cada banda 
ground_truth_array_raw = ds_ground_truth.GetRasterBand(1).ReadAsArray()
ground_truth_array = ground_truth_array_raw*filtered_mask
srs = ds_ground_truth.GetProjectionRef()
geo_transform = ds_ground_truth.GetGeoTransform()

cmap, norm = from_levels_and_colors([0,0.5,1,2],['white','#EA8A00','green'])

plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(111) ,plt.imshow(ground_truth_array, cmap=cmap),plt.title('Ground truth')
plt.show()

plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(111) ,plt.imshow(ground_truth_array[1400:1700,1200:1400], cmap=cmap),plt.title('Ground truth Zoom')
plt.show()

Ground truth details

In [18]:
def disease_info(array, geo_transform):
    unique, counts = np.unique(array[array!=0], return_counts=True)
    pix_width = geo_transform[1]
    pix_area = pix_width*pix_width

    lb_total_area = counts[0]*pix_area
    h_total_area = counts[1]*pix_area

    print("Late blight afected area: ", round(lb_total_area,2), "square meters")
    print("Healthy potato area: ", round(h_total_area,2), "square meters")
    print("Late blight vs Healthy potato area ratio: ", round(lb_total_area/h_total_area,2))
    print("Late blight area vs total area: ", round(lb_total_area/(lb_total_area+h_total_area),2))
    print("Healthy potato area vs total area: ", round(h_total_area/(lb_total_area+h_total_area),2))
In [19]:
disease_info(ground_truth_array, geo_transform)
Late blight afected area:  447.54 square meters
Healthy potato area:  348.48 square meters
Late blight vs Healthy potato area ratio:  1.28
Late blight area vs total area:  0.56
Healthy potato area vs total area:  0.44

Classification accuracy assessment

Random forest:

In [20]:
y_true = ground_truth_array.ravel()
target_names = ['0', 'Late blight', 'Healthy']

y_pred_rf = rf_classification.ravel()


print(classification_report(y_true, y_pred_rf, target_names=target_names))
              precision    recall  f1-score   support

           0       1.00      0.95      0.97   4214854
 Late blight       0.55      0.93      0.69    436912
     Healthy       0.87      0.67      0.76    340210

    accuracy                           0.93   4991976
   macro avg       0.81      0.85      0.81   4991976
weighted avg       0.95      0.93      0.93   4991976

Disease details of Random forest

In [21]:
disease_info(rf_classification, geo_transform)
Late blight afected area:  764.13 square meters
Healthy potato area:  267.9 square meters
Late blight vs Healthy potato area ratio:  2.85
Late blight area vs total area:  0.74
Healthy potato area vs total area:  0.26

Gradient Boost Classifier

In [22]:
y_pred_gbc = gbc_classification.ravel()

print(classification_report(y_true, y_pred_gbc, target_names=target_names))
              precision    recall  f1-score   support

           0       1.00      0.95      0.97   4214854
 Late blight       0.53      0.83      0.65    436912
     Healthy       0.71      0.67      0.69    340210

    accuracy                           0.92   4991976
   macro avg       0.75      0.82      0.77   4991976
weighted avg       0.94      0.92      0.92   4991976

Disease details Gradient Boost Classifier

In [23]:
disease_info(gbc_classification, geo_transform)
Late blight afected area:  705.11 square meters
Healthy potato area:  326.92 square meters
Late blight vs Healthy potato area ratio:  2.16
Late blight area vs total area:  0.68
Healthy potato area vs total area:  0.32

Support Vector Classifier

In [24]:
y_pred_svc = svc_classification.ravel()

print(classification_report(y_true, y_pred_svc, target_names=target_names))
              precision    recall  f1-score   support

           0       1.00      0.95      0.97   4214854
 Late blight       0.54      0.89      0.67    436912
     Healthy       0.83      0.70      0.76    340210

    accuracy                           0.92   4991976
   macro avg       0.79      0.84      0.80   4991976
weighted avg       0.95      0.92      0.93   4991976

Disease details Support Vector Classifier

In [25]:
disease_info(svc_classification, geo_transform)
Late blight afected area:  737.04 square meters
Healthy potato area:  294.99 square meters
Late blight vs Healthy potato area ratio:  2.5
Late blight area vs total area:  0.71
Healthy potato area vs total area:  0.29

Linear Support Vector Classifier

In [26]:
y_pred_lsvc = lsvc_classification.ravel()

print(classification_report(y_true, y_pred_lsvc, target_names=target_names))
              precision    recall  f1-score   support

           0       1.00      0.95      0.97   4214854
 Late blight       0.54      0.89      0.67    436912
     Healthy       0.82      0.70      0.76    340210

    accuracy                           0.92   4991976
   macro avg       0.79      0.85      0.80   4991976
weighted avg       0.95      0.92      0.93   4991976

Disease details Linear Support Vector Classifier

In [27]:
disease_info(lsvc_classification, geo_transform)
Late blight afected area:  733.6 square meters
Healthy potato area:  298.42 square meters
Late blight vs Healthy potato area ratio:  2.46
Late blight area vs total area:  0.71
Healthy potato area vs total area:  0.29

K-Neighbors Classifier

In [28]:
y_pred_knn = knn_classification.ravel()

print(classification_report(y_true, y_pred_knn, target_names=target_names))
              precision    recall  f1-score   support

           0       1.00      0.95      0.97   4214854
 Late blight       0.53      0.90      0.67    436912
     Healthy       0.82      0.65      0.72    340210

    accuracy                           0.92   4991976
   macro avg       0.79      0.83      0.79   4991976
weighted avg       0.95      0.92      0.93   4991976

Disease details K-Neighbors Classifier

In [29]:
disease_info(knn_classification, geo_transform)
Late blight afected area:  758.68 square meters
Healthy potato area:  273.35 square meters
Late blight vs Healthy potato area ratio:  2.78
Late blight area vs total area:  0.74
Healthy potato area vs total area:  0.26

References

Glasbey, C.A., 1993. An Analysis of Histogram-Based Thresholding Algorithms.

Otsu, N., 1979. A Threshold Selection Method from Gray-Level Histograms. IEEE transactions on systems, man, and cybernetics 9, 62–66.arXiv:1011.1669v3.

Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011